Structure formation during the collapse of a dipolar atomic Bose-Einstein condensate 
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We investigate the collapse of a trapped dipolar Bose-Einstein condensate. This is performed by 
numerical simulations of the Gross-Pitaevskii equation and the novel application of the Thomas- 
Fermi hydrodynamic equations to collapse. We observe regimes of both global collapse, where 
the system evolves to a highly elongated or flattened state depending on the sign of the dipolar 
interaction, and local collapse, which arises due to dynamically unstable phonon modes and leads 
to a periodic arrangement of density shells, disks or stripes. In the adiabatic regime, where ground 
states are followed, collapse can occur globally or locally, while in the non-adiabatic regime, where 
collapse is initiated suddenly, local collapse commonly occurs. We analyse the dependence on the 
dipolar interactions and trap geometry, the length and time scales for collapse, and relate our findings 
to recent experiments. 

PACS numbers: 03.75.Kk, 75.80.-|-q 



Wavepacket collapse is a phenomenon seen in diverse 
physical systems whose common feature is that they obey 
non-hnear wave equations e.g., in nonlinear optics 
0, plasmas 3 and trapped atomic Bose-Einstein con- 
densates (BECs) 0, H, E 0, S, In the latter case, 
collapse occurs when the atomic interactions are suffi- 
ciently attractive. For the usual case of isotropic s-wave 
interactions experiments have demonstrated both global 
and local collapse depending upon, respectively, 
whether the imaginary healing length is of similar size or 
much smaller than the BEG During global collapse 
the monopole mode becomes dynamically unstable and 
the BEG evolves towards a point singularity, with the 
threshold for collapse generally exhibiting a weak depen- 
dence on trap geometry [ll|, |T2|. Local collapse occurs 
when a phonon mode is dynamically unstable such that 
the collapse length scale is considerably smaller than the 
BEG. 

Recently, the Stuttgart group demonstrated collapse in 
a BEG with dipole-dipole interactions, where the atomic 
dipoles werepolarized in a common direction by an exter- 
nal field Q • The long-range nature of dipolar interac- 
tions means that the Gross-Pitaevskii wave equation that 
gov erns the BEG is not only non-linear but also non-local 
|l3l 0, [H, [l^. On top of being long-range, dipolar in- 
teractions are also anisotropic, being attractive in certain 
directions and repulsive in others. This anisotropy has 
manifested itself experimentally in the stability of the 
ground state, which is strongly dependent on the trap 
geometry [3], and in the anisotropic collapse of the con- 
densate Q . Some uncertainty exists over the mechanism 
of collapse in these systems. In the latter experiment, 
striking images indicate that the condensate underwent 
global collapse, which is li kely to have occurred through a 
quadrupole mode [l6l [TtI. IT8|. In the former experiment, 
however, recent theoretical results suggest that local col- 
lapse played a dominant role [Toj . 

A unique feature of trapped dipolar BEGs in compari- 



son to s-wave BECs is that they are predicted to exhibit 
minima in their excitation spectrum at finite values of 
the excitation quantum number [10, [HI, IH HI, [H, HI] , 
reminiscent of the roton minimum in the quantum liquid 
He-II. For gaseous BECs the depth of the minimum is 
tunable via microscopic parameters such as the dipolar 
and s-wave interaction strengths. An important physical 
consequence of the "roton" minimum is that it can lead 
to density modulations in the ground state dipolar BEG 
[IS [13, [3, [H [H, [13, [11]. However, the regions of pa- 
rameter space where they occur are small and lie close 
to the unstable region [Ij, [1^ . It is therefore natural to 
ask if these characteristic density modulations form dur- 
ing collapse where they might be more readily visible. 
As we shall see, collapse experiments can indeed provide 
us with an indirect yet accessible way of studying such 
effects. 

The dipolar BEGs that have been realized thus far 
[29! . [30j are formed of ^^Cr atoms with magnetic dipole 
moments d and coupling strength C'dd — fJ-ad^ , where 
fiQ is the permeability of free space. Ultracold quan- 
tum gases of polar molecules, which feature electric dipo- 
lar interactions, are also likely to be formed in the 
near future (sij . The dipolar interactions typically co- 
exist with s-wave interactions of characteristic amplitude 
g = ATrh^as/m, where as is the s-wave scattering length, 
and an important parameter is the ratio [3^ . 

Edd = Cdd/3g. (1) 

The s-wave coupling g can be effectively tuned between 
positive and negative infinity using a Feshbach resonance, 
and this was employed to control the stability of the dipo- 
lar BEG in the recent experiments [Sj] and Both the 
amplitude and sign of the dipolar coupling C'dd can also 
be tuned by rotation of the polarization axis [33| . A huge 
parameter space of dipolar interactions, from ejd = ~oo 
to -f 00, is thereby accessible for study. 

Insight into the collapse instability of a trapped dipolar 
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BEC can be obtained by considering its ground states. 
Such a route could be followed by beginning with a sta- 
ble ground state and adiabatically tuning the parame- 
ters towards the instability, which we term adiabatic col- 
lapse. Consider tuning the ratio £dd: as it is increased the 
ground state density profile evolves so as to benefit from 
the attractive part of the interaction. This can happen 
both globally and locally. For example, when Cdd > the 
system can undergo global magnetostriction and elongate 
along the axis of polarization [ll, [l^ IHj HI, Hll , tending 
towards a collapsed state of a line of end-to-end dipoles. 
However, the dipoles can also rearrange themselves lo- 
cally [13, m H m, m, [13, Hi and in particular, in a 
pancake-shaped geometry for Cdd > a "red blood cell" 
density profile has been predicted (23 |. 

If collapse is triggered suddenly, which we term non- 
adiabatic collapse, the system does not follow the ground 
state solutions and excitations play a role. In the dipo- 
lar collapse experiment j9i] the collapse was initiated by 
a change in g which took place over approximately one 
trap period. This time scale lies on the border between 
adiabatic and non-adiabatic collapse. 

Motivated by these issues we study theoretically the 
global and local collapse of a dipolar BEC over a sig- 
nificant range of edd that is accessible to current experi- 
ments. Our analysis is based on simulations of the Gross- 
Pitaevskii equation and the hydrodynamic (Thomas- 
Fermi) approximation, including the novel application 
of the hydrodynamic equations of motion to collapse. 
We observe that collapse occurs anisotropically with 
the dipoles tending to align themselves end-to-end for 
Cdd > and side- by-side for Cdd < 0. If collapse is 
approached adiabatically it occurs globally for moderate 
dipolar interactions (— l^edd^2) and beyond this we 
have also observed signs of local collapse. When collapse 
is initiated suddenly, it is dominated by the formation of 
local density structures, whose shape is determined by 
the dipolar interactions and can be related to unstable 
Bogoliubov modes. We map out the length and time 
scales for collapse, and the role of interaction strength 
and trap geometry. We then compare our results to re- 
cent experiments, in particular that of Lahaye et al. [^, 
where the condensate appeared to undergo global col- 
lapse. We show that our results are consistent with this 
observation and, furthermore, indicate how local collapse 
could be induced in this current experimental set-up. 

In Section U we introduce the mean-field description 
of dipolar BECs and, by employing the Thomas-Fermi 
approximation, derive the static solutions of the system 
and hydrodynamic equations of motion. In Section [Til we 
discuss the static solutions and their threshold for col- 
lapse. In Section IlIII we consider non-adiabatic collapse, 
induced by a sudden change in the interaction strength, 
and compare the hydrodynamic predictions with simula- 
tions of the Gross-Pitaevskii equation. In Section [TVl we 
extend our analysis of non-adiabatic collapse to map out 
the time and length scales of collapse, and in the latter 
case, show that the homogeneous Bogoliubov spectrum 



gives good agreement with the observed local collapse. In 
Section|V]we relate our findings to recent collapse experi- 
ments jSo^] and in Section lvTl we present our conclusions. 

I. THEORETICAL FRAMEWORK 

A. Dipolar Gross-Pitaevskii equation 

For weak interactions and at zero temperature the 
mean- field "wave function" of an atomic BEC tp = ip{r, t) 
satisfies the Gross-Pitaevskii equation (GPE), 

'^'^ = (-^^' + + 9\^\' + *dd(r, t)) V. (2) 
ot \ 2m J 

We assume that the confining potential has the 
cylindrically-symmetric harmonic form, 

V;xt(r) = imc^2(^2^^2^^2^2)^ (3) 

where ujx is the radial trap frequency and 7 = Wz/w^: is 
the trap's aspect ratio (which will be henceforth termed 
the trap ratio). Note that the trap has axial and ra- 
dial harmonic oscillator lengths Oz = y^h/mujz and Ox = 
y^h/mujx. Dipolar atomic interactions are described by 
the non-local mean-field potential <I>dd [H, [13, [H, [13 , 

*dd(r)= ydV[/dd(r-r')|^(r')p. (4) 

The interaction potential between two dipoles separated 
by r and aligned by an external field along a unit vector 
e is given by, 

C^dd (r) = — e,e, ^3 ■ (5) 

Throughout this work we consider the dipoles to be 
aligned in the ^-direction. It is useful to specify the dipo- 
lar interaction strength in Eq. ([2|) by the parameter, 

fcdd = iVfldd/a^, (6) 

where add = Cdd'7^/(127r/i^) is the dipolar "scattering 
length" and iV is the total number of condensate atoms. 
Note that the ratio £dd can be written as Edd = Odd/os- In 
the dipolar BEC collapse experiments 0, [^ k^A lies in the 
range 25 — 50. We will assume throughout this work that 
.g > such that the case of Edd < corresponds to Cdd < 
0. While the opposing case of 5 < is experimentally 
accessible, the s-wave interactions will typically induce 
collapse, rather than the dipolar interactions, and so will 
not be considered in this work. 

It is important to note that the basic GPE is insuffi- 
cient to describe the full collapse dynamics since higher- 
order effects, e.g. three-body loss, can become significant 
as the density escalates. However, the GPE provides an 
excellent prediction for the onset of collapse [1, [ll|, [l2| 
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and can be expected to accurately describe the early col- 
lapse dynamics. Our study will therefore consider the 
dynamics up to this point. Extension of these results to 
the full collapse dynamics could be made in future by 
including a three-body loss term in the GPE Ul] . 



B. Thomas- Fermi solutions 

Static solutions of the GPE can be expressed as 
'ip{r, t) = V'o(r) exp[— where ji is the chemical po- 
tential. We enter the Thomas-Fermi (TF) regime when 
the interaction energy dominates over the energy arising 
from density gradients, termed the zero-point kinetic en- 
ergy [4l[ . We then neglect the zero-point energy and the 
atomic density ng = IV'oP satisfies the expression, 



V{v) + gno{v) + ^dd{r) IJL 



(7) 



For an s-wave BEC, the ratio of interaction energy to 
zero-point energy is commonly specified as fcg — Nas/ux, 
with the system entering the TF regime when fcg 3> 1. 
The criterion for the TF regime is not so simple for a 
dipolar BEC since the anisotropic interactions make it 
strongly dependent on its shape. However, in the limit- 
ing cases of a highly elongated "cigar" condensate, or a 
highly flattened "pancake" condensate, the TF regime is 
valid when, respectively [36j . 



7^ 



dd 



fcdd 
73/2£dd 



(1 - £dd) > 1 [cigar] 



(1 -f 2£dd) > 1 [pancake]. 



(8) 
(9) 



A solution of Eq. ^ is given by an inverted parabola 
of the form 



and Rzo = Rxo/i^o- For an s-wave BEC, we retrieve the 
simple and expected result that Kq = 7- The anisotropic 
interactions, however, lead to magnetostriction such that 
kq < 7 for Edd > and kq > 7 for Sdd < 0. 

Consider the broader range of cylindrically-symmctric 
density profiles which have the same form as Eq. (fTO|) 
but are not limited to the equilibrium solutions. These 
general profiles define an energy "landscape" in terms of 
Rx and k given by [35| . 



E{Rx,k) 
N 



rriLulRl 
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15NgK 
28nRl 



1 - edd/(A^)](14) 



Then the static solutions pT|) correspond to stationary 
points in this landscape located at {Rxo, kq). By examin- 
ing the landscape in the vicinity of the solution one can 
determine whether it is stable (a global or local energy 
minimum) or unstable (maximum or saddle point). 

The TF regime can be formally regarded as the N — > 
00 limit of Gross-Pitaevskii theory. Indeed, the TF pre- 
dictions for the stability of a dipolar BEC depend only 
upon Edd and 7, thus simplifying the parameter space. 
However the TF model does not accurately describe sit- 
uations where the zero-point energy is considerable, e.g. 
close to the collapse threshold. In order to discern the 
effect of the zero-point energy we shall compare the TF 
solutions with those of a variational approach based upon 
a gaussian ansatz, as detailed in Appendix A, which in- 
cludes this energy contribution [l^, [3l • To describe de- 
viations from the TF (or gaussian) solutions, e.g. "red- 
blood cell" states [2?], one must solve the full GPE. This 
will also be considered in due course. 



C. Thomas-Fermi equations of motion for collapse 
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valid where no(r) > and no{r) = elsewhere, and the 
TF radii of the condensate are denoted by Rxo-, Ryo and 
Rzo- Note that the 0-subscript denotes static solutions. 
When the trap is cylindrically-symmetric about the same 
direction as the polarizing field, as we assume here, the 
density is also cylindrically symmetric. Its aspect ratio 
kq = Rxo I Rzo satisfies a transcendental equation [H,[33] , 



7^ 
where 



3edd/('«o) / 7 



- 2£dd - 1 



= edd-l. (11) 



To enable a hydrodynamic interp retation we employ 
the Madelung transform tp{r,t) = ^/n{r, t) exp[iS'(r, t)], 
where n{r,t) and S{r,t) are the density and phase dis- 
tributions, and define the "fluid" velocity as v{r,t) = 
{h/m)S/S{r,t) 38]. In the TF limit the dipolar GPE 
then leads to hydrodynamic equations (4lj ]. 



V + gn + $, 



dd 



(15) 
(16) 



Following earlier work for s-wave BECs [3^, Ho] there 
exists a class of exact time-dependent scaling solutions 
to Eqs. (HSl) and (HH) given by [H], 



/(«) 



1 + 3K^arctanh\/l — 



1-^2 (1 _ ^2)3/2 

The equilibrium radii are given by [33 ]. 
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Rl{t) Rl{t) Rl{t) 
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(17) 



-V [ax(t)x' + ay{t)y'' + az{t)z''\ (18) 



valid where n{r,t) > and n{r,t) = elsewhere, and 
no{t) = 15N/{8'KRxit)Ry{t)Rzit)) is the peak density 
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The time evolution of the radii Rj is governed by three 
ordinary differential equations, with the components of 
the velocity field given by aj = Rj/Rj, where j = x, y, z. 

Restricting ourselves to cylindrically-symmetric dy- 
namics where Rx{t) = Ry{t), and introducing the scaling 
factors Xx{t) = Rx{t)/Rxo and \z{t) — Rz{t)/Rzo (recall 
that the 0-subscript denotes the initial static solution), 
the time evolution is determined by the coupled ordinary 
differential equations. 



Xx — —ljXx + 



Xt.Xz 

1 

A2 



1 

3 2 fi'^0^x/Xz) 

2"o 4X1 ~ XI 



uh^Xz + 
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3 /(koXx/Xz) 
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where rj = 15N/47rmR^Q. These equations are signifi- 
cantly less demanding to solve than the full GPE and 
in certain limits analytic solutions exist (see Eq. ([26]) be- 
low). The TF equations of motion have been successfully 
applied to model the condensate dynamics under time- 
dependent trap ping including the important case of bal- 
listic expansion [4l|, \^ . They describe two independent 
collective excitation modes of the system: the monopole 
mode (when Xx{t) and Xz{t) are in phase) and the axis- 
symmetric quadrupole mode (when Xx{t) and Xz(t) are 
180° out of phase). In this paper we employ the TF 
equations of motion to study global collapse. In the pure 
s-wave case global collapse occurs through the monopole 
mode [l3| , but when dipolar interactions dominate global 
collapse occurs through the quadrupole mode [igI. [TtI. [lH . 

We will employ the TF equations of motion to describe 
non-adiabatic collapse, triggered by a sudden change in 
SAd{t)- To conform to the current experimental method 
P, 9] we shall implement this through a sudden change 
in the s-wave interactions g(t). By, (i) starting with the 
BEG well below the threshold for collapse and (ii) sud- 
denly changing to a state which is well above the thresh- 
old for collapse, the regime where zero-point energy domi- 
nates can be bypassed and the TF equations should apply 
throughout. We will verify this in due course. 



II. ADIABATIC COLLAPSE 

Imagine an experiment that starts with a stable ground 
state BEG and adiabatically increasing the magnitude of 
£dd- The condensate will remain in the ground state cor- 
responding to the instantaneous value of ejd and will 
eventually collapse at a critical value of e^d [H, [iBl • The 
threshold for collapse in general depends upon £dd, 9, 
7 and N. We consider two possible scenarios for adi- 
abatic collapse: i) if there is no roton minimum in the 
excitation spectrum then adiabatic collapse proceeds in a 



similar manner to the usual pure s-wave case, i.e. a global 
collapse via a low-lying shape oscillation mode once the 
(imaginary) healing length becomes of the same order as 
the condensate size [lo| : ii) if there is a roton minimum 
then this deepens as £dd increases and, at the point at 
which it reaches zero energy, can lead to local collapse 
on a length scale determined by the roton minimum. 
We note that technically speaking neither scenario can 
be truly adiabatic since the mode responsible for col- 
lapse has zero frequency at the collapse threshold, but 
the non-adiabaticity can be confined to a small region of 
parameter space if Edd is increased slowly enough. 



A. Global adiabatic collapse 

Ground state solutions, characterised by their aspect 
ratio kqj are presented in Fig.[T]Ja) as a function of Edd for 
various trap ratios 7. The figure shows the predictions 
given by numerical solution of the GPE, the parabolic 
TF solution, and the gaussian variational ansatz. We 
first consider the parabolic TF solutions (solid lines in 
Fig. [TJa)), which, we recall, can be characterized solely 
by Edd and 7. For Edd = we observe that kq =7, as 
expected for s-wave condensates. We will consider the 
regimes of positive and negative Sdd separately with the 
aid of typical energy landscapes pictured in Fig. [T]^b): 
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FIG. 1; (Color online) (a) Aspect ratio no of the ground state 
solutions as a function of Sdd for various trap ratios 7. Pre- 
sented are the stable TF solutions of Eq. (solid lines), 
GPE solutions with kdd = 115 (crosses), and variational so- 
lutions for fcdd = 115 (dashed line), (b) Energy landscapes of 
Eq. (Ull) for 7 = 2 and £dd = (i) -0.8, (ii) -0.52, (iii) 0.8, (iv) 
1.025 and (v) 1.4. Light/dark regions correspond to high/low 
energy, and white contours help to visualise the landscapes. 
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• Edd > 0: As Edd is increased from zero, kq de- 
creases since the dipoles prefer to lie end-to-end 
along z. The experiment [4J] observed this magne- 
tostrictive effect, in good agreement with the TF 
predictions. For < Edd < 1, the parabolic solu- 
tions (llOj) are global minima of the TF energy func- 
tional ^ (Fig. [T](b)(iii)). For Sdd > 1, however, 
the TF solution becomes only a local minimum 
(Fig. [Hb)(iv)), with the global minimum being a 
collapsed state of zero width. Indeed, there is an 
upper critical dipolar-to-s-wave ratio, e'^ > 1, be- 
yond which the local minimum disappears and the 
whole system is unstable to collapse into a kq = 
state (Fig. [ljb)(v)), i.e. an infinitely thin line of 
dipoles. Note that e[j^ depends on the trap ratio 7. 
Elongated BECs, being predominantly attractive, 
are most unstable, with e'^ ~ 1. Increasing 7 in- 
creases due to the increasing repulsive interac- 
tions in the system. Indeed, for 7 > 5.17, e^J = 00 
[l^[ll, UHl, i.e., the parabolic solutions are robust 
to collapse for any interaction strength (see the case 
of 7 = 6 in Fig. [Ha)). 

• Edd < 0: As £dd is decreased from zero, kq in- 
creases since the dipoles now prefer to lie side-by- 
side in the transverse plane. For —0.5 < Edd < 0, 
the parabolic solutions are robust to collapse, while 
for £dd < —0.5, they become local energy minima 
with the global minima being a collapsed plane of 
dipoles (Fig. [Hb)(ii)). We define a critical value 

< —0.5, below which no stable parabolic solu- 
tions exist (Fig.[T](b)(i)) and the system is unstable 
to collapse into a kq = 00 state, i.e., an infinitely 
thin plane of side-by-sidc dipoles. Pancake geome- 
tries are particularly prone to this with e|^j « —0.5, 
while the system becomes increasingly stable in 
more elongated geometries. Indeed, in sufficiently 
elongated geometries 7^0.19, the parabolic solu- 
tions are stable to collapse with ej^^ — —00. 

Numerical solutions of the time-independent GPE for 
kdd — 115 are shown in Fig. [Ija) as crosses. Our method 
of determining the BEG widths is detailed in Appendix 
B. While the GPE solutions are generally in good agree- 
ment with the TF results, deviations become significant 
near the point of collapse where the zero-point energy 
has a considerable stabilising effect on the solutions. We 
have additionally performed time-dependent simulations 
in which Edd is varied slowly, and find that the conden- 
sate follows these ground state solutions and maintains 
a parabolic-like density profile. This is consistent with 
global collapse. We have performed a similar analysis for 
kdd = 14 and observe the same qualitative results, with 
the threshold for collapse pushed to slightly higher values 
of l^ddl due to enhanced zero-point energy. 

The dashed lines in Fig. [TJa) are the results given by 
the gaussian variational method of Appendix A. Note 
that if zero-point effects were neglected, the gaussian so- 
lutions would satisfy the same transcendental equation 



dTl]) as the parabola, i.e., the aspect ratio of the BEG k 
is largely independent of the ansatz (although the pro- 
file and radii do differ). Indeed, for small values of |edd| 
the variational and TF predictions are almost identical. 
However, close to the onset of collapse these predictions 
deviate. Importantly, the gaussian variational method 
gives excellent agreement with the full GPE solutions 
right up to the point of collapse, highlighting the impor- 
tance of zero-point energy at the point of collapse. 



B. Local adiabatic collapse 

A surprising result of the parabolic density profiles is 
that within the TF approximation they remain stable 
even as Sdd — > 00 and £dd —00, providing the trap 
ratio is sufficiently extreme (7 > 5.17 for edd > and 
7 < 0.19 for Edd < 0). Similar behaviour arises for a 
variational approach based on a gaussian density pro- 
file [1^, [1^ . The recent experiment of Koch et al. Q has 
observed the stability of a dipolar condensate under vari- 
ous trap ratios, in reasonable qualitative agreement with 
the gaussian and TF predictions. However, numerical 
solutions of the GPE show that, even though the stabil- 
ity is enhanced under extreme trap ratios, collapse will 
occur for strong enough dipolar interactions [1^|. This 
apparent contradiction arises because the gaussian and 
parabolic solutions are only capable of describing low- 
lying monopole and quadrupole fluctuations. In the pres- 
ence of a roton minimum the instability of a mode with 
high quantum number can lead to local adiabatic col- 
lapse. Bohn et al. recently employed a theoretical model 
that allowed for local collapse and found improved quan- 
titative agreement with the experimental observations in 
pancake geometries [l^, suggesting that local collapse 
does indeed play a key role in such geometries. Note 
that, by contrast, in s-wave condensates an adiabatic re- 
duction in g should always induce global collapse. 

In order to induce local adiabatic collapse it is neces- 
sary to use values of £dd that fall outside — 1 ^ £dd ^ 2 
(the region of Fig. 1). We have probed pancake-shaped 
ground state solutions that extend beyond this range. 
For 7 = 10 and 20 we have probed up to Sdd = 30 with 
no evidence for collapse. However, for 7 = 6 we have ob- 
served the adiabatic onset of a local collapse instability 
at Edd ~ 2.2, characterised by the formation and growth 
of cylindrical density shells (similar to those that will be 
discussed in Section 11111) . Our results are consistent with 
those of Ronen et al. [24|], who predicted that, in a purely 
dipolar BEG (edd = 00), collapse is possible for large 7 
and that, close to the collapse threshold, the ground state 
adopts density corrugations. It is reasonable to assume 
that as one passes into the unstable region that these lead 
to local collapse. This picture was also recently suggested 
by Bohn et al. [l^. Note that for the moderate range 
of £dd that we concentrate on here, the TF solutions are 
a very good approximation to the Gross-Pitaevskii solu- 
tions and the adiabatic collapse proceeds globally. 
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III. NON-ADIABATIC COLLAPSE 

An alternative way to induce collapse is to perform a 
non-adiabatic change of of the form, 

We will assume that this is achieved by tuning the s-wave 
interactions from go — Cdd/Se^^j to gf = Cdd/Se^d- 

In order for the TF equations of motion to remain valid 
we require that the TF approximation holds throughout 
the dynamics i.e. regimes where the zero-point energy 
is important are avoided. This requires that the initial 
condensate is well inside the regime of stable solutions 
(IejjjI < l^ddD' ^^'^ subsequently that the system is deep 
within the collapse regime (|eddl > l^ddD- 
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FIG. 2: (Color online) Collapse dynamics in a cigar trap 7 = 
0.1 for = 0.8 and = 1.4. (a) Energy landscape of 
Eq. (fT4|l at t = 0. (b) Energy landscape for t > 0, with the 
ensuing TF trajectory indicated (green/grey line). Energy 
is scaled in units of hcux- (c) Xx(t) and (d) Xz{t) from the 
full TF equations of motion (black solid line), the simplified 
case of Eqs. (|23I24|) (solid grey line), and CPE simulations for 
kdd = 80 (grey dotted line) and 2000 (black dotted line), (e) 
Evolution of the energy components during CPE simulations 
with fcdd = 80 (grey lines) and k^d = 2000 (black lines). 
Shown are the kinetic (solid line), zero-point E^p (dotted 
line) and dipolar Ed (dashed line) energies, renormalised by 
the total energy Etot- 



A. Non-adiabatic collapse for e^d > 

We begin with a stable condensate with e[|j^ = 0.8 and 
suddenly switch to e^^ > 1. 



1. Cigar trap 

For a cigar trap 7 = 0.1 the initial state lies in the 
energy landscape of Fig. [2ta). For t > 0, we switch to 
= 1.4 and an energy landscape of Fig. [DJb). This 
system is unstable and the BEC undergoes a trajec- 
tory (green/grey line) towards the collapse {Rx = 0) 
region. According to the TF scaling parameters Xx{t) 
and Xz{t) (solid black hues in Fig. O^c) and (d)), the 
condensate accelerates to zero width in the x-direction 
after t « 0.75u;~^, during which the width reduces by a 
very small amount, less than 1%. The collapse is there- 
fore highly anisotropic. Note that this collapse occurs 
relatively fast since the condensate is initially in an elon- 
gated state, close to the collapse threshold. As we will see 
the time for collapse is strongly dependent on the initial 
shape of the condensate and therefore on the trap ratio, 
with more pancake condensates taking longer to collapse. 

In our example the BEC is highly elongated both ini- 
tially and throughout its dynamics. Assuming K{t) <C 1 
analytic results for the TF equations of motion can then 
be obtained. Expanding /(k) as. 



- = 1 + 4k2 + log(«;/2) + ©(k^), (22) 



(a) (i) (ii) (iii) 




-4 -2 2 4 -4 -2 2 4 -4 -2 2 4 
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FIG. 3: (Color online) Snapshots of the collapse of the cigar 
BEC with fedd = 2000 in the (a) x-y plane and (b) y-z plane 
at (i) t = 0, (ii) t = 0.47 and (iii) t = 0.57a;~\ 
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then Eqs. (|T9|) and ([20)) become, to lowest order, 



~ -ujIXx + [1 - £dd(i)] 
\, ~ -t^^T^A^ + [1 - £dd(i)] 



A3 A, 



(23) 
(24) 



Recall that 77 = 15iV/47rr7ii?i^Q. Analogous equations have 
been derived for an expanding repulsive s-wave BEC [40| . 
For our time-dependent protocol (pij) and to lowest order 
in kq, Eq. (p4)) . has the solution, 



A2(t) = cos(cj3;7t). 



(25) 



This corresponds to the limit of the non-interacting gas, 
and the next correction is of order Kg. To zeroth order 
in 7, Xz{t) = 1, and Eq. has the solution. 



-j=\/il + cr) cos{2uJxt) + l-a (26) 
v2 



where 



(27) 



These simplified analytic solutions can give a remark- 
ably good description of the dynamics. For example, in 
Fig. IH^c), their prediction of Xx{t) is in excellent agree- 
ment with the full TF equations of motion. For Xz{t) 
(Fig. m^d)) deviations are clearly visible, although the 
dynamics are very slow in this direction. 

We have performed GPE simulations of the collapse for 
a BEC with fcdd — 80. These results (grey dotted lines 
in Fig.[2Ic) and (d)) are not in good agreement with the 
TF predictions. This is because the TF condition for 
an elongated dipolar BEC ^ is not remotely satisfied 
[3^ . While the BEC approximates a TF profile in the 
z-direction the transverse profile is more akin to a non- 
interacting gaussian ground state. Under a much larger 
interaction strength of fcdd = 2000 (black dotted line), 
which does satisfy the TF criterion ([5]), we find good 
agreement with the TF predictions up until t w 0.65a;~^. 
We have evaluated the energy contributions to the GPE 
as outlined in Appendix B and plotted them in Fig. [2je) 
for both fcdd = 80 (grey lines) and 2000 (black lines). 
The validity of the TF approximation for fcdd = 2000 is 
confirmed by the smallness of the zero-point energy in 
comparison to the other energy contributions. Indeed, it 
remains small up until t w O.IlUx, thereby validating the 
use of the TF approach. At t = the dipolar energy 
-Ed is negative, indicating the attractive configuration of 
dipoles in the initial state. The dipolar energy remains 
negative and grows in magnitude as collapse proceeds. 
Meanwhile the total kinetic energy grows and diverges 
at i « 0.707^:^. The point at which the energies diverge 
effectively marks the breakdown of the validity of the 
numerical simulations and the TF approach. 

In Fig. [3] we present snapshots from the GPE simula- 
tions with fcdd — 2000. The initial density [Fig. EKi)] is 



highly elongated along z and approximates the TF in- 
verted parabola. Up to t ~ 0.4a;~^ the BEC collapses 
anisotropically while maintaining the inverted parabola 
shape. However, from this point in time [Fig. El^ii)] a lo- 
cal density structure emerges (highlighted by contours), 
characterised by modulations in the density. This struc- 
ture evolves to form a striking arrangement of ellipsoidal 
shells of high and low density [Fig.[31Jiii)], with a high den- 
sity region in the centre of the condensate. The regions 
of high density grow in population and peak density, and 
thereby undergo collapse. Note that the development of 
the density wave structure marks a clear deviation of the 
system from the parabolic TF solutions. 

According to the Bogoliubov spectrum of a homoge- 
neous system (to be discussed further in Section IIV[) 
the dipolar BEC can be dynamically unstable to phonon 
modes [l3|, [35| . This instability has a strong dependence 
on angle relative to the polarization direction, and phys- 
ically this represents the tendency for the system to pref- 
erentially align the dipoles in an end-to-end configuration 
in the polarization direction. This effect gives rise to the 
highly elongated shell structures observed here. 



2. Pancake trap 

We now consider a BEC within a pancake trap 7 = 2. 
The initial state with — 0.8 resides in the energy 
landscape of Fig. [H^a). Following the sudden switch to 
= 1.4 and an energy landscape of Fig. lUJb), the BEC 
follows a trajectory (green/grey line) towards the collapse 
Rx = region. According to the TF equations of motion 
(solid lines in Fig.HJJc) and (d)), the condensate acceler- 
ates to zero width in the x-direction after t « 2a;~^, while 
the z-width oscillates by approximately 10%. Again, the 
collapse is highly anisotropic. It is considerably slower 
than in the cigar trap because the initial condensate is in 
a more stable state, dominated by repulsive interactions. 
Note that if the condensate begins in a highly fiattened 
state K ^ 1, it elongates over time and its aspect ratio 
will be reversed. As such no expansion of k, analogous 
to Eqs. (^51) and ([M)) . is appropriate. 

For fcdd = 80 the TF pancake criterion ^ is satis- 
fied. The corresponding GPE predictions (dotted lines 
in Fig. djc) and (d)) agree well with the TF predictions 
up to i « 2a;~^, with the zero-point energy remaining 
small up until this point. Initially, the dipolar energy 
is positive due to the dominance of repulsive dipoles in 
the initial fiattened BEC but as the BEC collapses, it 
becomes negative and diverges. Again, collapse causes a 
divergence in the kinetic energy. 

During the collapse we again see the formation of a 
local density structure, as shown in Fig. [5fa,b). These 
structures are not closed ellipsoidal shells but now recti- 
linear shells aligned in the z-direction. In this case it is 
likely that the large trap frequency in z, which leads to 
a large excitation energy for axial excitations, suppresses 
the curvature of the shells. Here the structure features a 
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FIG. 5: (Color online) Snapshots of the BEC density during 
the collapse of Fig. |3]in (a) x-y plane and (b) y-z plane at (i) 
t = 0, (ii) t = IloZ^ and (iii) t = l.SuJ^. Light/dark regions 
correspond to high/low density, while contours in (ii) help to 
visualise the density structures. 



anti-node, rather than a node. 



FIG. 4: (Color online) Collapse dynamics in a pancake trap 
7 = 2 for = 0.8 and e^d = (a) Energy landscape of 
Eq. ifT^ at t = 0. (b) Energy landscape for t > 0, with the 
TF trajectory indicated (green/grey line). Energy is scaled in 
units of hiox- (c) \x(t) from the TF equations of motion (solid 
line) and the GPE with fcdd = 80 (dotted line), (d) Same for 
Az(t). (e) Kinetic Ey, zero-point _Ezp and dipolar Ed energies 
from the GPE, renormalised by the total energy i?tot- 



B. Non-adiabatic collapse for e'l^ < 



Here we start with a stable condensate with eJI^ 



and suddenly switch to < — 0.5 



-0.2 



node of low density for r = and a single high density 
shell at finite radius. This structure evolves rapidly with 
atoms moving away from the centre of the condensate to 
populate the outer shell. This causes the observed flat- 
tening in \x at t « 2cl'^^. By later times [Fig. (SJiii)] the 
shell has elongated axially and shrunk radially, and the 
system continues to undergo local collapse. 

To probe the effect of the trap geometry we also present 
the collapse dynamics for 7 = 5. For e^^ = 1.4 no col- 
lapse occurs and so we employ a more extreme value of 
= 4. The TF equations of motion (solid lines in 
Fig- Eta) and (b)) predict that the condensate undergoes 
shape oscillations rather than collapse. The GPE results 
(dotted lines) agree up to t = l.huj~^. During this time 
the condensate density (Fig. [S]Jd-e)(i)) remains approxi- 
mately an inverted parabola, although weak local density 
perturbations begin to emerge (Fig.[6l^d-e)(ii)). However, 
for t > l.bu!^^, the kinetic energy i?k diverges to -l-oo and 
diverges to —00. We again observe the development 
(and collapse) of local density structure (Fig.[6Ud-e)(iii)). 
This structure features considerably more shells than our 
previous example, primarily due to the larger radial ex- 
tent of the condensate. It also features a central density 



In an elongated trap 7 = 0.2 we have observed collapse 
for e^^ — —4 with the dynamics presented in Fig. [T] Un- 
der the TF equations of motion (solid lines) , the conden- 
sate undergoes large oscillations in the x-direction and 
collapses slowly in z. The collapse occurs predominantly 
in the z direction towards an infinitely thin plane of side- 
by-side dipoles. This is the opposite to the regime of 
edd > 0. Note that the collapse is slow because the con- 
densate is initially in an elongated state where the attrac- 
tive interactions that induce collapse are weak. With 
fcdd = 80 the GPE initial state is in the TF pancake 
regime ([5]). The GPE results agree well up to their point 
of validity at i « StJ^^. During this time, the zero-point 
kinetic energy remains relatively small, and the kinetic 
and dipolar energies undergo large oscillations due to the 
radial shape oscillations. At t w Scj^^, i?k and di- 
verges signifying the limit of validity of the simulations. 

Consideration of the condensate density profile during 
collapse reveals that the condensate develops weak planar 
density corrugations by i « 4w~^ [Fig. [7{d,e)(ii)]. These 
become amplified into a striking planar density pertur- 
bation of approximately 10 localised planes of dipoles, 
aligned in the x — y plane [Fig. [7I^d,e)(iii)]. This is the 
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FIG. 6: (Color online) Dynamics in a pancake trap 7 = 5 
under a sudden change from edd ~ 0.8 to 4. (a) Xx{t) and 
(b) \z{t) from the TF equations (solid line) and the GPE 
with fcdd = 80 (dotted line), (c) Kinetic E]^, zero-point E^p 
and dipolar Ed energies from the GPE simulations, (d)-(e) 
Density snapshots during the GPE simulations in the (d) x-y 
plane and (e) y-z plane at (i) t = 0, (ii) 1 and (iii) 1.5lo^^ . 



same phenomenon as observed earlier but with perpen- 
dicular orientation due to the fact that the £dd < 
dipoles are now attractive when side-by-side. 



2. Pancake trap 

For = —0.8, the condensate dynamics within a pan- 
cake trap 7 = 2 are presented in Fig. [5] According to the 
TF equations of motion, A^, (t) accelerates to zero within 
t Q.6uj~^ while Xx{t) decreases more slowly. Taking 
^dd — 80 the initial ground state satisfies the pancake 
TF criteria (O and we see excellent agreement with the 
TF predictions up to t « 0.6w~^, during which Ezp re- 



mains relatively small. However, beyond this point E^p 
diverges, as does and E^- 

Like the cigar case, the condensate density develop pla- 
nar density perturbations, as shown in Fig.[8l^d,e)(ii) and 
(iii). However, given that the pancake system is narrower 
in z, only two planes of high density emerge. 



IV. GENERAL PROPERTIES OF COLLAPSE 



We now map out the general behaviour of the impor- 
tant length scales and time scales of the collapse, and the 
critical trap ratio that can stabilise against collapse. 
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FIG. 7: (Color online) Collapse dynamics in a cigar trap 7 = 
0.2 for e% = -0.2 and Edd = -4. (a) A^(t) and (b) A^t) 
from the TF equations of motion (solid line) and the GPE 
with fcdd = 80 (dotted line).(c) Kinetic i5k, zero-point E^p 
and dipolar E^ energies, (d)-(e) GPE density in the (d) x-y 
plane and (e) y-z plane at (i) t = 0, (ii) 4.4 and (iii) 5ujx^ ■ 
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FIG. 8: (Color online) Collapse dynamics in a pancake trap 
7 = 2 under a sudden change from Sdd = —0.2 to —0.8. 
(a) Xx{t) and (b) Xz{t) according to the full TF equations 
of motion (solid line) and CPE simulations with kdd = 80 
(dotted line).(c) Kinetic i?k, zero-point E^p and dipolar Ed 
energies from the CPE, renormalised by the total energy Etot ■ 
(d)-(e) Density snapshots in the (d) x-y and (e) y-z planes at 
(i) t = 0, (ii) 0.42 and (iii) 0.53lu~\ 



A. Length scales for collapse 

A homogenous dipolar BEC is unstable to periodic 
density perturbations (phonons) when Sdd > 1 or Edd < 
—0.5 [ij]. This can be seen immediately from the Bo- 
goliubov dispersion relation between the energy Eb and 
momentum p for phonons in the gas, given by, 



^--V(y +2,n{l + ...(3cos^^-l)}|^. 

(28) 

A mode with energy evolves as exp{iEBt/H) and so 
when Eb becomes imaginary, the mode grows exponen- 



tially, i.e., a dynamical instability. This dispersion rela- 
tion depends on 0, the angle between the momentum 
of the phonon and the external polarizing field. For 
Edd > 1, phonons propagating perpendicularly to the po- 
larization direction {9 — 7r/2) undergo a dynamical insta- 
bility while no instability occurs along the polarization 
direction {6 = 0). This is illustrated in Fig. [Hb) (in- 
set) which plots the spectrum for Sdd = 1-4 and for 
the polarizations 9 — (dashed line) and 6 = Tr/2 (solid 
line); the point where the solid line touches zero marks 
the transition to instability. In an infinite and initially 
homogeneous system we expect this instability to break 
the condensate up into a lattice of filaments, i.e. cylin- 
drical structures aligned along the polarization direction. 
In trapped condensates, shells or planes (aligned along 
z) may be favored instead. Meanwhile, when edd < —0.5 
(we remind the reader that we mean in this case that 
Cdd < 0), it is phonons propagating along the polar- 
ization direction (9 = 0) that undergo a dynamical in- 
stability while no instability occurs in the perpendicular 
direction {9 = 7r/2). This suggests an instability towards 
stratification into planes lying perpendicular to the po- 
larization axis. These predictions are consistent with our 
observations in the previous section. 

We can estimate the characteristic length scale of the 
collapse structures from the homogeneous Bogoliubov 
spectrum (|28p . If Pc is the critical momentum for which 
the dispersion relation passes through zero energy (which 
signifies the onset of dynamical instability), then we can 
expect the characteristic length scale to be Ic = 2-Kh/pc- 
Although we are concerned with trapped, inhomogeneous 
condensates, the same length scale Ic will apply providing 
the system is of size R ^ Ic- From Eq. it is straight 
forward to show that this length scale is given by. 



Ic = A 
L ^ B 



V^dd - 1 

_VM_ 

V|l + 2£dd| 



for edd>l (29) 
for Edd < -0.5. (30) 



Here £, — 1/ \/87rno|as| is the s-wave healing length at the 
condensate centre. A and B are factors that take account 
of the trapping, and are unity in a homogenous system. 

We can improve the applicability of our length scale 
predictions to inhomogeneous systems by taking into ac- 
count the trapping in one or two directions. In an in- 
finite cylindrical BEC with tight radial trapping, such 
that the radial density profile is a gaussian, it is known 
that the axial speed of sound is reduced by a factor of 
\/2 in comparison to a uniform system with the same 
peak density ng , 43:]. This rescaling of the dispersion re- 
lation arises because the average density is in fact given 
by no/2. In general, we can obtain the average density 
for both cigar and pancake condensates by integrating 
out the tightly confined direction. We will perform this 
for the relevant cases of gaussian and parabolic (TF) den- 
sity profiles. The resulting modification to the dispersion 
relations along the weakly confined direction(s) gives the 
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FIG. 9: Length scale of collapse Ic as a function of e^^j ac- 
cording to GPE simulations (circles) and the predictions of 
Eq. (|31|l and (I32|) for a homogeneous system A — B = 1 
(dashed lines) and for a trapped system A B 1 (solid 
line), (a) The regime Sdd < —0.5, assuming a cigar-shaped 
BEG with 7 = 0.2, e% = -0.2 and fcdd = -80. For the cigar 
prediction (solid line) we assume a TF transverse profile for 
which B — \/3/2. (b) The regime £dd > 1, assuming a pan- 
cake geometry 7 = 5 with e^d = 0.8 and fcdd = 80. For the 
pancake prediction (solid line) we assume a TF axial profile 
for which A = 1^/5/4. Inset: Homogeneous Bogoluibov spec- 
trum of Eq. (l28|l for Edd = 1.4 and for S = (dotted line) and 
9 — n/2 (dot-dashed line). Note that the error bars in the 
GPE results arise from the grid discretization. 



trapping parameters A and B which occur in Eq. ((29|) 
and ([50)1 . For edd > 1 the interesting case is that of a 
pancake, for which A Apau ss = 2^/"* w 1.2 in the gaus- 
sian case and A —i- Atf = '\/5/4 w 1.1 in the TF case. 
For £dd < —0.5 the interesting case is that of a cigar, for 
which B i?Gauss = \/2 w 1.4 in the gaussian case and 
B Btf = \/372 ~ 1-2 in the TF case. 

In order to relate our numerical results on adiabatic 
and non-adiabatic collapse to the unstable modes of 
Eq. psp. let us consider the case of positive £dd- For 
£dd = 1 the collapse length scale is infinite. If Edd is in- 
creased adiabatically then Ic decreases. At some critical 
£dd, Ic becomes equal to the condensate size R and a dy- 
namical instability occurs on the length scale of R. Thus, 
according to Eq. (p8)) . if the point of collapse is reached 
adiabatically, global collapse will occur. This picture is 
qualitatively consistent with our observations in Section 
Ullin the range — 1 < edd ^ 2, where we observed that adi- 
abatic collapse proceeds in a global manner. However, 
outside of this regime we saw evidence for adiabatic lo- 
cal collapse, which we attributed to the presence of a 
roton minimum in the excitation spectrum. The homo- 
geneous dispersion relation (p8|) does not have a roton 
minimum and in order to introduce one it is necessary 
to explicitly include a trap in at least one dimension. In 
highly cigar-shaped or pancake-shaped systems the ro- 
ton minimum occurs in the dispersion relation for low- 
energy excitations along the weakly confined direction. 
The momentum at which the roton minimum occurs, Pr, 
defines a length scale Ir = 2'Kh/pr which is usually closely 



related the system size in the tightly trapped direction 
(23I [23 |. Adiabatic local collapse is therefore expected 
to take place on length scales Ir where the roton mini- 
mum touches the zero energy axis. We have compared 
simple analytic predictions for the dispersion relations of 
an infinite cigar [22| and an infinite pancake [2^, which 
do include a roton minimum, with the homogeneous re- 
sult pS)) . We find that the value of Ic can significantly 
differ from the predictions of (|29)) and ((30)) when the sys- 
tem is close to the collapse threshold, but quickly become 
very similar as we move deeper into the collapse regime. 
Thus, we expect that an adiabatic collapse experiment 
could provide clear evidence for the presence of a roton 
minimum in the excitation spectrum, but a non-adiabatic 
collapse experiment would be less conclusive. To tackle 
this problem properly one should numerically obtain the 
full Bogoliubov spectrum of a trapped system (37| . 

We now consider the situation where we suddenly 
switch from < 1 to e^^j ^ 1. Here we go from the 
regime where Ic ^ R to Ic < R, and a dynamical in- 
stability occurs on a local scale. Although the dynamical 
instability evolves for all length scales in the range I > Ic, 
the imaginary energy eigenvalue is largest for / = Ic and 
this unstable mode dominates the system. This is consis- 
tent with our observations in Section lTlTl where, following 
a sudden change in Cdd, collapse evolves mainly through 
local density structures. 

We now specifically consider the case where collapse is 
induced suddenly by modifying the s-wave interactions, 
while the dipolar interactions remains constant. Then we 
can rewrite Eqs. ((29| and ((30|) as. 

These predictions are plotted in Fig. [D^a) and (b), 
for the regimes of edd < —0.5 and edd > 1, respec- 
tively. We have conducted a series of GPE simulations 
for trapped dipolar BECs to determine the true collapse 
length scale (defined as the distance between peaks in ad- 
jacent shells), with these predictions being shown by the 
circles. Note that in Fig. [DJa) we employ a cigar-shaped 
BEG with 7 = 0.2 and e^J^ = -0.2, and in Fig. Mh) we 
employ a pancake-shaped BEG with 7 = 5 and e°j = 0.8. 
Even in the homogeneous limit A — B — 1 (dashed lines) 
the analytic predictions for Ic are in very good agreement 
with the simulations. Examination of the GPE solutions 
reveals that the density profile in the tightly confined di- 
rection is closely approximated by the TF profile. The 
analytic res ults using the a ppro priate trapping parame- 
ters (A = ■\/5/4 and B = •\/3/2) are shown by solid lines 
in Fig. [H With these trapping parameters included the 
agreement becomes excellent, and clearly demonstrates 
the importance of taking the trapping into account. 
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B. Time for global collapse 

We now make some simple predictions for the charac- 
teristic collapse time Tc- Having demonstrated that the 
TF model gives a good description of collapse we will 
employ it exclusively here. Since the TF equations of 
motion cannot describe local collapse, our analysis is lim- 
ited to global collapse. Furthermore we will only consider 
non-adiabatic collapse since the time scale for adiabatic 
collapse is technically infinite. 

Experimentally, the collapse time is the time at which a 
'sudden' depletion of the condensate occurs due to three- 
body loss [9| . The suddenness arises because of the dra- 
matic scaling of losses with time: the rate of three-body 
loss scales as n? , and the density n itself accelerates in 
amplitude during collapse. Consequently, a good esti- 
mate for Tc is the time over which the peak density di- 
verges or, equivale ntly , the time over which one or more 
radii tend to zero UM ,- 

The TF collapse time can be obtained, in general, by 
numerical solution of the TF equations of motion. How- 
ever, for highly elongated BECs, Eqs. ((25)) leads to an 
analytic form for Tc given by. 



1 

- arccos 
2 



a - 1 



(33) 



where a is given by Eq. (I??)) . It should be noted that in 
oo, Eq. (I27p becomes. 



the limit 



dd 



lim 



_ Wo^o 

— ^dd 



(34) 



where we recall that 77 = / AirmR^Q. Thus the limit- 
ing value of Tc is determined only by the initial parame- 
ters 7 and £jj. 

In Fig. [TOTa') and (b) we show how Tc depends on 
the final interaction parameter e^^ for initial values of 
= —0.2 and 0.8, respectively. For weak interactions 
Tc diverges as the interactions become too weak to induce 
collapse, while in the limit of large interactions (positive 
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FIG. 10: Collapse time Tc following a sudden change in Sdd 
according to the TF equations of motion. Various trap ratios 
are presented, (a) e^d = —0.2 and edd < 0. (b) ejd — 0.8 
and £dd > 0. For the case of 7 = 0.1 we also plot the analytic 
expression of Eq. l|33)l (grey dashed line). 
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FIG. 11: Critical trap ratio 7c that stabiUses the BEC against 
collapse as a function of initial interaction strength e^di ac- 
cording to the full TF equations of motion. The final inter- 
action strength ejjd = 10* is so large that it is effectively the 
infinite limit. 



or negative edd), Tc tends towards a finite value, as ex- 
pected. For Edd < 0, elongated systems become more 
stable and are the slowest to collapse. In contrast, for 
Edd > 0, elongated systems are least stable and therefore 
the fastest to collapse. Note that for the highly elongated 
case of 7 = 0.1, the analytic collapse time of Eq. ((33)) is in 
excellent agreement with the full TF equations of motion. 



C. Critical trap ratio for global collapse 

In the range of £dd considered in this work, there exists 
a critical trap ratio 7c which can stabilise the parabolic 
TF solutions against global collapse. As we showed in 
Section HI, if collapse is approached adiabatically this 
critical trap ratio is fixed, being jc — 5.2 for positive edd 
and 7c = 0.19 for negative edd- However, if collapse is 
induced non-adiabatically, the critical trap ratio becomes 
a function of e^^^. To examine this threshold for global 
collapse in more detail. Fig. [TTj plots the critical trap 
ratio 7c as a function of the initial interaction parameter 



ejjd for both negative and positive edd- Note that we 
have considered je^^jl = 10^, which is so large that it 
effectively behaves like the infinite limit. For e°j < 0, 
7c tends towards zero for e°j — > 0, and in the opposing 
limit of edd — *■ — 00 it increases asymptotically towards 
the adiabatic value 7c = 0.19. Conversely, for e^j^ > 0, 
7c diverges as e^^ —> 0, and in the opposing limit of 
^dd ^ 00 it decreases towards the static value 7c = 5.17. 
Note that, providing |eddl ^ 0, 7c varies only weakly 
with eJljj and becomes very close to the static value of 7c. 



V. DISCUSSION 

In a recent series of experiments the Stuttgart group 
demonstrated the collapse of a dipolar condensate [sj, 
[^. A Feshbach resonance was employed to give time- 
dependent control over Os and therefore edd- In Ref. [1], 
as was reduced slowly (over several trap periods) to probe 
the stability of the ground state to collapse. For 7 = 1, 



13 



(a) (i) (ii) (iii) 



« 



mm 


5 


mm 


5 


















-5 




-5 




□ y o 

X (units of a^) 




^ ^ 

X (units of aj 




R 
U 

X (units of a J 


mm 


5 




6 




□ 














-5 




-5 




-5 5 
y (units of a,) 




-5 5 
y (units of aj 




-5 5 
y (units of aj 



FIG. 12: (Color online) Collapse dynamics in Stuttgart sys- 
tem under the linear decrease of as over time Tr — 1ms. Den- 
sity in the (a) x-y plane and (b) y-z plane at (i) t = 0, (ii) 
t = 3.1 and (iii) t = 3.7uj-^ . 



collapse was observed for Sdd ~ 1, which is in good agree- 
ment with the static solutions in Fig. [Ija). For 7 « 10, 
the BEC was stabilised under purely dipolar interactions, 
i.e. the limit Sdd 00, in good qualitative agreement 
with our understanding of the role of trap geometry. 

In Ref . Q , collapse was observed in an almost spherical 
trap. From an initial BEC at t = with 0° = 1.59nm, 
the scattering length was reduced linearly to a final value 
ttg = 0.27nm over a time scale r ~ 1ms. The BEC was ob- 
served to collapse anisotropically towards a narrow cigar 
shape, aligned in z, over a time scale of t°^p w 1.5ms. 
We will see that this time scale is sufficiently long that 
local collapse is not induced. Following this initial col- 
lapse an explosion occurred, resulting in a spectacular 
state with a shape resembling that of a d-wave orbital. 



Numerical simulations of the GPE with three-body loss 
were in excellent agreement with the observed dynamics 
i. 

The experimental collapse appears to occur globally 
with no signs of local collapse. To resolve the issue of 
the apparent absence of local collapse, we have simu- 
lated the experimental dynamics based on iV = 20, 000 
atoms and a fully anisotropic trap {uJx,ujy,u!z) — 2tt x 
(660, 400, 530)Hz. Under linear ramping of ag over time 
T = 1ms, our results in Fig. [12] confirm that the con- 
densate collapses globally rather than locally. The con- 
densate collapses, mainly in the x — y plane, and forms 
a very narrow cigar-shaped BEC. Because lu^ > i-^y, the 
transverse collapse is quickest in the x-direction. The 
elongated collapsed state (Fig. [121 iii)) forms after a col- 
lapse time Tc_~ 1ms. This is quicker than the time scales 
reported in [9[ , where it is known that eddy currents sup- 
press the applied change in as and effectively extend the 
ramping time by a factor of two or three. This does 
not affect our qualitative results but merely slows down 
the collapse. Indeed, if we employ a ramping time of 
T = 2ms, we find Tc ~ 1.6ms, in agreement with [9]. 

The experimental ramping time is far from being sud- 
den since it is of the order of the trap period that charac- 
terises the internal dynamics of the BEC. We can there- 
fore expect that global collapse will be initiated well be- 
fore any local instabilities. Indeed, if we make a more 
sudden change of as, e.g., t = 0, then we see in Fig. [T3| 
that local collapse now occurs. However, within this non- 
cylindrically symmetric geometry we observe the forma- 
tion of parallel density stripes, rather than shells. Again, 
the x-direction collapses towards zero width. However, 
the y-direction does not shrink globally but develops a 
corrugated structure that enables the dipoles to predom- 
inantly align along z. These stripes become amplified 
and collapse themselves. 
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FIG. 13: (Color online) Collapse dynamics in Stuttgart sys- 
tem under the sudden decrease in as. Density in the (a) x-y 
plane and (b) y-z plane at (i) t — 0, (ii) t — 1.08 and (iii) 
t = 1.16tJj\ 



VI. CONCLUSIONS 

We have studied the collapse dynamics of a dipolar 
Bose-Einstein condensate, triggered either by an adia- 
batic or non-adiabatic change in the dipolar-to-s-wave 
ratio Edd — Cdd/35. In general, the collapse occurs 
anisotropically and is driven by the dipoles seeking to line 
up end-to-end for edd > and side- by-side for Sdd < 0. 
In the case of adiabatic collapse, where the ground state 
solutions are followed up until the instant of collapse, 
we observe both global and local collapse. In the range 
— 1 ^ £dd ^ 2 we find global collapse towards a single line 
or plane of dipoles. Outside of this range we have seen 
adiabatic local collapse which we suggest is a signature 
of a roton minimum in the excitation spectrum. Similar 
theoretical predictions have also recently been reported 
by Bohn et al. [l^. Note that care must be taken to 
distinguish such local collapse from the results of non- 
adiabatic collapse. 

If collapse is triggered non-adiabatically via a sudden 
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change in edd, the instabihty can jump to length scales 
much less than the condensate size, resulting in local col- 
lapse. We have analysed this instability over the con- 
siderable range —10^ edd ^20. This instability can be 
understood in terms of the amplification of dynamically 
unstable phonon modes. For a cylindrically-symmetric 
condensate, the system develops a periodic structure of 
density shells for Edd > or disks for Sdd < 0, which 
become amplified and subsequently collapse. 

We applied the TF equations of motion, previously 
used to model oscillatory and expansion dynamics, to 
study non-adiabatic collapse of the condensate. We 
showed that this method is valid providing that zero- 
point kinetic effects remain small throughout. This can 
be ensured by employing a TF condensate initially, and 
suddenly switching the interactions deep into the collapse 
regime. The predictions are in excellent agreement with 
the full Gross-Pitaevskii equation for the majority of the 
collapse dynamics. The TF predictions fail when the 
condensate develops local density structures and thereby 
deviates from the TF parabolic density profile. 

Our results are consistent with the experiment by La- 
haye et al. @ . There the increase of Edd was sufficiently 
slow to ensure that global collapse dominates the system. 
However, by changing the interactions more suddenly, it 
should be possible to induce local collapse of the conden- 
sate. A rich variety of structures can be formed depend- 
ing upon the orientation of the dipoles and the aspect 
ratio of the trap. These structures are related to those 
that have been predicted to occur in the ground state, 
but during local collapse they become amplified. Such 
transient structures could be observed experimentally by 
Bragg scattering of light. 



APPENDIX A: GAUSSIAN VARIATIONAL 
SOLUTIONS 

When the TF approximation is not valid (e.g. for weak 
interactions) , it can be appropriate to employ a gaussian 
ansatz. We consider a cy lindrically-symmetric gaussian 
ansatz of the form [1, [l^ . 



I kN 



exp 



(Al) 



where is the transverse width and k is the aspect ratio. 
The energy of this ansatz is. 
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Variational solutions are then obtained by minimising 
this with respect to Ix and k. 

APPENDIX B: NUMERICAL SOLUTION OF 
THE DIPOLAR GPE 

A split-step fast fourier transform method is employed 
to evolve V'lr) on a three-dimensional spatial grid [37|, 



typically with 64'^ points. The dipolar atomic potential 
in A;-space is given by. 



f>dd(fc) 



4^Cdd 



(Bl) 



where kx^ ky and kz are the cartesian wavevectors and 
k'^ = kx + ky + k'^. By combining Eq. ^ with the convo- 
lution theorem, the dipolar mean-field potential is given 
by $dd(r) = T^^[Udd{k)n{k)], where the inverse Fourier 
transform, denoted by J^~^ , is performed numerically us- 
ing a standard fast Fourier transform algorithm. For pan- 
cake and spherical systems we employ the corrections to 
Udd{k) outhned in Ref. [sj. Ground-state solutions of 
the dipolar GPE are found by propagation in imaginary 
time t —it. From a suitable initial guess and with 
renormalisation at each time step, the GPE converges to 
the ground state of the system. 

To monitor the size of the BEG we evaluate (z^) and 
(x^) and relate them to the TF radii Rz and Rx via. 



(z2) = J drz^m' = ^Rl 
(x^) = J dr.'\M' - jRl 



(B2) 



To compare the numerical widths directly with the TF 
scaling parameters Xz{t) and Xx{t) we define, 



(B3) 
(B4) 



Furthermore, the total energy i^tot of the dipolar BEG 
is numerically evaluated from the GPE energy functional 
[4H and given by. 




Etot - (l^l^^l' + ^-tl^l' + *dd 



(B5) 

The terms in the integrand correspond to, from left to 
right, the kinetic energy i?k, potential energy Ep, dipo- 
lar interaction energy Ed and s-wave interaction energy 
Eg ■ Note that the zero-point kinetic energy is the kinetic 
energy associated with variations in the real part of ip. 
Expressing ■0 = V'r + where ip^ and ipi are real quan- 
tities, the jVV'p-term in Eq. (|B5p can be decomposed as 
jVi/'p = [V'f/'r]^ + [VV'i]2. The zero-point kinetic energy 
can thus be defined as. 



E., 



2m 



(B6) 
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